﻿using MathNet.Symbolics;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

using MathNet.Numerics.LinearAlgebra;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.Statistics;
//using MathNet.Symbolics;
using Expr = MathNet.Symbolics.SymbolicExpression;
using MathNet.Numerics.Distributions;
using System.Diagnostics;
using System.Net;
using System.Xml.Linq;
using static System.Runtime.InteropServices.JavaScript.JSType;
using System.Xml;
using MathNet.Numerics;

namespace NumericalAnalysis.DataInterpolation
{
    public class NewtonDifferenceInterpolation
    {

        public Expr polynomial = 0;   //多项式表达式
        List<double> polyCoefficients = new List<double>(); //各项系数 幂次高至低
        List<int> coefficienOrders = new List<int>(); //各项幂次

        static Vector<double> theVec = Vector<double>.Build.Dense(2);
        List<Vector<double>> vectors = new List<Vector<double>>();


        public Matrix<double> diffInterpMat { get; set; }

        List<double> xs = new List<double>();
        List<double> ys = new List<double>();


        public string interpType { get; set; }

        int count = 0;

        /// <summary>
        /// 初始化
        /// </summary>
        /// <param name="_vectors"></param>
        public NewtonDifferenceInterpolation(IEnumerable<Vector<double>> _vectors, string _interpType = "forward")
        {
            this.interpType = _interpType;

            count = _vectors.Count();
            vectors = _vectors.ToList();

            xs = _vectors.Select(v => v[0]).ToList();
            ys = _vectors.Select(v => v[1]).ToList();

        }

        public NewtonDifferenceInterpolation(IEnumerable<Tuple<double, double>> _tuples, string _interpType = "forward")
        {
            this.interpType = _interpType;

            count = _tuples.Count();
            vectors = _tuples.Select(v =>
            {
                var vec = Vector<double>.Build.Dense(2);
                vec[0] = v.Item1;
                vec[1] = v.Item2;
                return vec;
            }).ToList();

            xs = _tuples.Select(v => v.Item1).ToList();
            ys = _tuples.Select(v => v.Item2).ToList();


        }


        public Matrix<double> DiffInterpolationForward()
        {
            diffInterpMat = Matrix<double>.Build.Dense(count, count, 0.0d);
            //零阶差分
            for (int i = 0; i < count; i++)
            {
                diffInterpMat[i, 0] = ys[i];
            }

            for (int j = 1; j < count; j++)
            {
                for (int i = 0; i < count - j; i++)
                {
                    diffInterpMat[i, j] = diffInterpMat[i + 1, j - 1] - diffInterpMat[i, j - 1];
                }
            }
            return diffInterpMat;
        }

        public Matrix<double> DiffInterpolationBackward()
        {
            diffInterpMat = Matrix<double>.Build.Dense(count, count, 0.0d);
            //零阶差分
            for (int i = 0; i < count; i++)
            {
                diffInterpMat[i, 0] = ys[i];
            }

            for (int j = 1; j < count; j++)
            {
                for (int i = j; i < count; i++)
                {
                    diffInterpMat[i, j] = diffInterpMat[i, j - 1] - diffInterpMat[i - 1, j - 1];
                }
            }
            return diffInterpMat;
        }

        /// <summary>
        /// 常规公式解法
        /// </summary>
        /// <param name="interpType">forward   backward</param>
        /// <returns></returns>
        public Expr FitInterpolation()
        {
            var t = Expr.Variable("t");
            //polynomial = 0;
            if (polynomial.ToString().Contains("t")) return polynomial;

            if (interpType != "forward" && interpType != "backward") return null;
            if (xs.IsEquidistance() == 0) return null;

            List<double> dCoeff = new List<double>();
            Matrix<double> diffInterp = Matrix<double>.Build.Dense(count, count, 0);
            if (interpType == "forward")
            {
                diffInterp = DiffInterpolationForward();
                dCoeff = diffInterp.Row(0).ToList();

                double factorial = 1;
                Expr termPoly = t;  //差分项
                polynomial = dCoeff[0];
                for (int i = 1; i < count; i++)
                {
                    polynomial += dCoeff[i] * termPoly / factorial;
                    termPoly *= (t - i);
                    factorial *= (i + 1);
                }
            }
            if (interpType == "backward")
            {
                diffInterp = DiffInterpolationBackward();
                dCoeff = diffInterp.Row(count - 1).ToList();

                double factorial = 1;
                Expr termPoly = t;  //差分项
                polynomial = dCoeff[0];
                for (int i = 1; i < count; i++)
                {
                    polynomial += dCoeff[i] * termPoly / factorial;
                    termPoly *= (t + i);
                    factorial *= (i + 1);
                }

            }

            //显示多项式
            Console.WriteLine("--------------------------------");
            Console.WriteLine("多项式 polynomial:");
            Console.WriteLine(polynomial);

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式简化 polynomial-ExponentialSimplify:");
            //Console.WriteLine(polynomial.ExponentialSimplify());

            //Console.WriteLine("--------------------------------");
            //Console.WriteLine("多项式展开式 polynomial-Expand:");
            //Console.WriteLine(polynomial.Expand());

            polynomial = polynomial.Expand();

            return polynomial;
        }


        /// <summary>
        /// 从给出的向量计算
        /// </summary>
        /// <param name="sourceVecs"></param>
        /// <returns></returns>
        public IList<Vector<double>> CalculateBy(IList<Vector<double>> sourceVecs)
        {
            for (int i = 0; i < sourceVecs.Count(); i++)
            {
                double xValue = sourceVecs[i][0];
                var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
                double yValue = polynomial.Evaluate(symbols).RealValue;
                //double   yValue = CalculateBy(xValue);
                sourceVecs[i][1] = yValue;
            }
            return sourceVecs;
        }

        /// <summary>
        /// 从给出的数列计算
        /// </summary>
        /// <param name="xValues"></param>
        /// <returns></returns>
        public IList<double> CalculateBy(IList<double> xValues)
        {
            var symbols = new Dictionary<string, IEnumerable<double>>
                { { "t", xValues } };
            List<double> yValues = polynomial.CalculateBy(symbols).ToList();
            return yValues;
        }

        /// <summary>
        /// 从给出的数计算
        /// </summary>
        /// <param name="xValue"></param>
        /// <returns></returns>
        private double CalculateBy(double xValue)
        {
            //复杂多项式避免使用Compile方法计算  将严重拖慢性能
            //Func<double, double> func = polynomial.Compile("t"); 
            //double yValue = func(xValue);

            var symbols = new Dictionary<string, FloatingPoint> { { "t", xValue } };
            double yValue = polynomial.Evaluate(symbols).RealValue;

            return yValue;
        }

        /// <summary>
        /// 计算
        /// </summary>
        /// <param name="xTerm">终端x值  前向时为首个 ，后向为末个</param>
        /// <param name="xValues">要求点的X</param>
        /// <param name="dis">插值点的等距步长</param>
        /// <returns></returns>
        public IList<double> CalculateBy(double xTerm, IList<double> xValues, double dis)
        {
            var symbols = new Dictionary<string, IEnumerable<double>>
                { { "t", xValues.Select(x=>(x-xTerm)/dis) } };
            List<double> yValues = polynomial.CalculateBy(symbols).ToList();

            return yValues;
        }

        public IList<double> CalculateT(List<double> interpXs, List<double> checkXs, string type)
        {
            var xTerm = 0d;
            if (type == "forward") xTerm = checkXs.First();
            if (type == "backward") xTerm = checkXs.Last();
            double xDis = interpXs[1] - interpXs[0];

            var symbols = new Dictionary<string, IEnumerable<double>>
                { { "t", checkXs.Select(x=>(x-xTerm)/xDis) } };
            var checkYs = polynomial.CalculateBy(symbols);
            return checkYs.ToList();
        }


    }


    public class NewtonDifferenceInterpolationTestDemo
    {
        string var_t = "t";
         
        //public IList<double> CalculateT(List<double> interpXs, List<double> checkXs, string type)
        //{


        //    var xTerm = 0d;
        //    if (type == "forward") xTerm = checkXs.First();
        //    if (type == "backward") xTerm = checkXs.Last();
        //    double xDis = interpXs[1]  - interpXs[0] ;

        //    var symbols = new Dictionary<string, IEnumerable<double>>
        //        { { "t", checkXs.Select(x=>(x-xTerm)/xDis) } };
        //    var checkTs = polynomial.CalculateBy(symbols);


        //    var checkdPnts = new List<Tuple<double, double>>();
        //    for (int i = 0; i < checkXs.Count; i++)
        //    {
        //        var tuple = new Tuple<double, double>(checkXs[i], checkTs[i]);
        //        checkdPnts.Add(tuple);
        //    }


        //    return checkTs.ToList();
        //}

        /// <summary>
        /// sin函数
        /// </summary>
        /// <returns></returns>
        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleSin(string interpType = "forward")
        {
            List<double> interpXs = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.1);
            List<double> checkXs = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.02);   //构造sin函数插值点 

            var t = Expr.Variable(var_t);
            Expr sourceExpr = t.Sin();  //原函数


            //原函数的离散点
            var interpPnts = sourceExpr.GeneratePointBy(var_t, interpXs).ToList();

            ////进行插值  
            var interp = new NewtonDifferenceInterpolation(interpPnts, interpType);
            var expression = interp.FitInterpolation();  //插值函数
            var checkYs = interp.CalculateT(interpXs, checkXs, interpType);

            var checkdPnts = new List<Tuple<double, double>>();
            for (int i = 0; i < checkXs.Count; i++)
            {
                var tuple = new Tuple<double, double>(checkXs[i], checkYs[i]);
                checkdPnts.Add(tuple);
            }

            return (interpPnts, checkdPnts);
        }

        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleRunge(string interpType = "backward")
        {
            List<double> interpXs = NumericalService.LinSpace(-5, 5, 0.5);
            List<double> checkXs = NumericalService.LinSpace(-4, 4, 0.4);

            var t = Expr.Variable(var_t);
            Expr sourceExpr = "1 / (t * t + 1)";
             
            //原函数的离散点
            var interpPnts = sourceExpr.GeneratePointBy(var_t, interpXs).ToList();

            ////进行插值  
            var interp = new NewtonDifferenceInterpolation(interpPnts, interpType);
            var expression = interp.FitInterpolation();  //插值函数
            var checkYs = interp.CalculateT(interpXs, checkXs, interpType);

            var checkdPnts = new List<Tuple<double, double>>();
            for (int i = 0; i < checkXs.Count; i++)
            {
                var tuple = new Tuple<double, double>(checkXs[i], checkYs[i]);
                checkdPnts.Add(tuple);
            }

            return (interpPnts, checkdPnts);

        }

        public (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleHomeWork(string interpType = "forward")
        { 
            List<Tuple<double, double>> interpPnts = new List<Tuple<double, double>>();
            List<double> values = new List<double>()
                { 12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 };
            for (int i = 0; i < values.Count; i++)
            {
                double xk = i * 2;
                double yk = values[i];
                Tuple<double, double> tuple = new Tuple<double, double>(xk, yk);
                interpPnts.Add(tuple);
            }
            var interpXs = interpPnts.Select(v => v.Item1).ToList();
            var interpYs = interpPnts.Select(v => v.Item2).ToList();

            List<double> checkXs = NumericalService.LinSpace(0 + 0.1, 24 - 0.1, 0.1);
             
            var interp = new NewtonDifferenceInterpolation(interpPnts, interpType);

            var interpExpr = interp.FitInterpolation();

            var checkYs = interp.CalculateT(interpXs, checkXs, interpType);
            var checkdPnts = new List<Tuple<double, double>>();
            for (int i = 0; i < checkXs.Count; i++)
            {
                var tuple = new Tuple<double, double>(checkXs[i], checkYs[i]);
                checkdPnts.Add(tuple);
            } 
            return (interpPnts.ToList(), checkdPnts.ToList());
             
        }

    }




    //public class NewtonDifferenceInterpolationTest
    //{
    //    /// <summary>
    //    /// 测试
    //    /// </summary>
    //    /// <param name="exsistPnts">存在插值点</param>
    //    /// <param name="needPntXs">需要求的点自变量</param>
    //    /// <returns></returns>
    //    public static IList<Tuple<double, double>> NewtonDifferenceInterpolationCalc(IList<Tuple<double, double>> exsistPnts, IList<double> needPntXs)
    //    {
    //        var t = Expr.Variable("t");
    //        NewtonDifferenceInterpolation interp =
    //            new NewtonDifferenceInterpolation(exsistPnts, "backward");

    //        var expression = interp.FitInterpolation();

    //        //var needPntYs = interp.CalculateBy(needPntXs);

    //        var xTerm = 0d;
    //        if (interp.interpType == "forward") xTerm = needPntXs.First();
    //        if (interp.interpType == "backward") xTerm = needPntXs.Last();
    //        double xDis = exsistPnts[1].Item1 - exsistPnts[0].Item1;

    //        var needPntYs = interp.CalculateBy(xTerm, needPntXs, xDis);

    //        var tuples = new List<Tuple<double, double>>();
    //        for (int i = 0; i < needPntXs.Count; i++)
    //        {
    //            var tuple = new Tuple<double, double>(needPntXs[i], needPntYs[i]);
    //            tuples.Add(tuple);
    //        }
    //        return tuples;
    //    }


    //    public static (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleHomeWork()
    //    {
    //        var t = Expr.Variable("t");
    //        List<Tuple<double, double>> extVecs = new List<Tuple<double, double>>();

    //        List<double> values = new List<double>()
    //            { 12, 9, 9, 10, 18, 24, 28, 27, 25, 20, 18, 15, 13 };
    //        for (int i = 0; i < values.Count; i++)
    //        {
    //            double xk = i * 2;
    //            double yk = values[i];
    //            Tuple<double, double> tuple = new Tuple<double, double>(xk, yk);
    //            extVecs.Add(tuple);
    //        }
    //        var interpYs = NewtonDifferenceInterpolationCalc(extVecs, extVecs.Select(v => v.Item1).ToList()).ToList();

    //        List<double> array2 = NumericalService.LinSpace(0 + 0.1, 24 - 0.1, 0.1);
    //        var checkdYs = NewtonDifferenceInterpolationCalc(extVecs, array2).ToList();

    //        return (interpYs, checkdYs);
    //    }


    //    /// <summary>
    //    /// sin函数
    //    /// </summary>
    //    /// <returns></returns>
    //    public static (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleSin()
    //    {
    //        var t = Expr.Variable("t");
    //        List<Tuple<double, double>> extVecs = new List<Tuple<double, double>>();
    //        //构造sin函数插值点
    //        List<double> xs = NumericalService.LinSpace(0, Math.PI * 4, Math.PI * 2 * 0.1);
    //        extVecs = NumericalService.CalcByFunc((x) => { return Math.Sin(x); }, xs).ToList();

    //        ////进行插值 并 返回检查点    这里插值点就是检查点  
    //        var interpYs = NewtonDifferenceInterpolationCalc(extVecs, extVecs.Select(v => v.Item1).ToList()).ToList();
    //        interpYs = extVecs;
    //        List<double> array2 = NumericalService.LinSpace(0, Math.PI * 2, Math.PI * 2 * 0.02);
    //        var checkdYs = NewtonDifferenceInterpolationCalc(extVecs, array2).ToList();

    //        return (interpYs, checkdYs);
    //    }


    //    /// <summary>
    //    /// 龙格函数
    //    /// </summary>
    //    /// <returns></returns>
    //    public static (List<Tuple<double, double>>, List<Tuple<double, double>>) ExampleRunge()
    //    {
    //        var t = Expr.Variable("t");
    //        List<Tuple<double, double>> extVecs = new List<Tuple<double, double>>();
    //        //构造runge函数插值点
    //        Func<double, double> rungeFunc = (x) => { return 1 / (x * x + 1); };
    //        List<double> array = NumericalService.LinSpace(-5, 5, 0.5);
    //        extVecs = NumericalService.CalcByFunc(rungeFunc, array).ToList();
    //        ////进行插值 
    //        var interpYs = NewtonDifferenceInterpolationCalc(extVecs, extVecs.Select(v => v.Item1).ToList()).ToList();

    //        List<double> array2 = NumericalService.LinSpace(-4, 4, 0.02);
    //        var checkdYs = NewtonDifferenceInterpolationCalc(extVecs, array2).ToList();

    //        return (interpYs, checkdYs);
    //    }

    //}

}
